Cite this Article
Tang Yao-Zong, Li Xiao-Lin. Meshless analysis of an improved element-free Galerkin method for linear and nonlinear elliptic problems . Chinese Physics B, 2017, 26(3): 030203  
Permissions
Meshless analysis of an improved element-free Galerkin method for linear and nonlinear elliptic problems
Tang Yao-Zong1, Li Xiao-Lin2, †
College of Mathematics and Statistics, Kashgar University, Kashgar 844000, China
College of Mathematics Science, Chongqing Normal University, Chongqing 400047, China

 

† Corresponding author. E-mail: lxlmath@163.com

Abstract
Abstract

We first give a stabilized improved moving least squares (IMLS) approximation, which has better computational stability and precision than the IMLS approximation. Then, analysis of the improved element-free Galerkin method is provided theoretically for both linear and nonlinear elliptic boundary value problems. Finally, numerical examples are given to verify the theoretical analysis.

PACS: 02.60.Cb;02.60.Lj;02.30.Em
Keyword:meshless method;moving least squares approximation;element-free Galerkin method;error estimate
1. Introduction

Over the past half century, the numerical solutions of many physics and engineering problems have been obtained by mesh-based methods such as the finite element method. In these methods, the approximation of unknown variables is related exactly to the geometry of the elements. Meshless (or meshfree) methods,[1,2] in which the approximation of unknown variables requires only nodes, can overcome the meshing-related drawbacks.[1–3]

Shape functions play an important role in meshless methods. The moving least squares (MLS) approximation is one of the most extensively used methods to construct meshless shape functions.[4–6] In recent years, some MLS variants,[1] such as the complex variable MLS,[7] the interpolating MLS,[4,8–10] the complex variable interpolating MLS,[11,12] and the improved interpolating MLS,[13–15] have also been presented. The MLS approximation is developed from the traditional least squares method (LSM), which means that the LSM is used at each computing point in practical numerical processes. A disadvantage of the LSM is that the final system of algebra equations is sometimes ill-conditioned. Hence, the MLS may form an ill-conditioned or singular moment matrix. Besides, the matrix must be inverted, which leads to the decrease in stability and the increase in computational cost and error. To overcome this drawback, weighted orthogonal basis functions are generated for the MLS.[16,17] Afterwards, the improved moving least squares (IMLS) approximation is developed. A key advantage of the IMLS approximation is that the moment matrix is diagonalized and then the system of algebra equations is solved without deriving the inverse matrix. Consequently, the burden of inverting the moment matrix is eliminated, and the IMLS approximation has a higher computing efficiency.[17–19]

Recently, a shifted and scaled monomial basis function has been developed to stabilize the MLS.[20] Theoretical analysis and numerical verification show that both the determinant and the condition number of the moment matrix in the stabilized MLS are independent of the nodal distance, and thus the stabilized MLS prevents the instability occurrence.[21] The inherent instability of the interpolating MLS has also been studied theoretically and numerically.[22] We should point out that the shifted and scaled monomial basis function has already been widely used in another class of meshless method, i.e., the reproducing kernel particle method (RKPM), since 1995.[23,24] Over the past two decades, there have been many developments on the RKPM. The equivalence between the MLS and the kernel approximations has been addressed.[25,26] Actually, when the monomial basis function is used, the element-free Galerkin (EFG) method and the RKPM are equivalent.

In this paper, shifted and scaled orthogonal basis functions are developed to stabilize the IMLS approximation. The determinant and the condition number of the moment matrix in the stabilized IMLS are invariable with respect to the nodal distance. Therefore, the stabilized IMLS is expected to obtain more accurate and stable results.

The improved element-free Galerkin (IEFG) method is a typical IMLS-based meshless method.[1] In this method, the IMLS is used to approximate the unknown variables, and the Galerkin weak form is employed to obtain the discretized system equations. The IEFG method has been successfully applied to problems in fracture,[27] elastic,[28] transient heat conduction,[29] elastodynamics,[30] biological population,[31] and viscoelasticity,[32] and to wave equation,[33] Camassa and Holm equation,[34] and Schrödinger equation.[35] It is shown that, under the same node distribution, the numerical error of the IEFG method is much less than that of the EFG method; and under the same accuracy, the IEFG method needs fewer nodes than the EFG method. Therefore, the IEFG method has greater computational precision and efficiency than the EFG method.

The IEFG method has gained great success in computational physics and engineering. Contrarily, only a few published papers deal with the associated mathematical theory. In Ref. [36], Li et al. analyzed theoretically the properties and errors of the IMLS, and then deduced the error estimates of the IEFG method for linear Poisson equations. In Refs. [37] and [38], Li et al. established error estimates for the MLS, and then obtained error analysis of the Galerkin boundary node method. For the EFG method, error estimates have been established for the linear Poisson equation,[1,21,39,40] elasticity equation,[41] sine-Gordon and sinh-Gordon equations.[42]

The second goal of this paper is to analyze theoretically the IEFG method for both linear and nonlinear elliptic equations. Detailed computational formulas are deduced. Besides, based on the error results of the IMLS, the theoretical error analysis of the IEFG method for linear and nonlinear elliptic equations is provided in energy norms.

The rest of this paper is outlined as follows. Section 2 presents and analyzes a stabilized IMLS approximation. Then, analysis of the IEFG method for both linear and nonlinear elliptic problems is presented in Sections 3–6. Finally, numerical examples and conclusions are given in Sections 7 and 8, respectively.

2. A stablized IMLS approximation
2.1. Formulations

Let be a set of N nodes in a bounded domain (n = 1,2,…). As in the MLS,[4] the local approximation ℐu of a function u in the IMLS can be defined as

(1)
where is an approximation operator, are constructed basis functions, are unknown coefficients, m is the number of the basis, and is the evaluation point x or a nodal point .

The basis used in the MLS is a complete set of linear independent monomial functions. In the IMLS, the basis is composed of a complete set of weighted orthogonal functions that can be constructed by the Schmidt orthogonalization technique as

(2)
where the inner product about nodes is defined for functions and as[17]
where , , are nodes with influence domains that cover the point x, denotes the global sequence numbers of these nodes, and are weight functions.

In Eq. (2), the given function can be chosen as the monomial basis functions used in the MLS.[16] Let , then in this case a linear basis is given by

and a quadratic basis by
Besides, we can define as[17]
where for one-dimensional problems, and or for n-dimensional problems.

The coefficients in Eq. (1) are obtained by minimizing the norm

as
where with are the nodal values, and

Since the functions given by Eq. (2) are orthogonal, i.e., for i ≠ j, we have

(3)

Substituting Eq. (3) into Eq. (2) and setting , the global approximation of u(x) is

(8)
where the IMLS shape function is

2.2. Choice of basis functions

As in Refs. [21] and [22], it can be proved theoretically that the determinant and the 2-norm condition number of the moment matrix A(x) in the IMLS are

(4)
(5)
where
is the nodal spacing, and are constants independent of h, and the notation denotes the largest degree of the functions .

Equations (4) and (5) indicate that, for h small enough, A(x) in the IMLS becomes poorly conditioned and the accuracy of computing decreases.

In Refs. [21] and [22], a shifted and scaled monomial basis function is used in the MLS and the interpolating MLS to improve the condition number of the associated moment matrix. In this study, the following shifted and scaled functions are defined using as

where is a fixed point located on the influence domain of the evaluation point x. In particular, we can take . Then, as in Eq. (2), shifted and scaled orthogonal basis functions are constructed using as !

Finally, substituting for in the IMLS, a stabilized IMLS can be deduced directly. The corresponding moment matrix is denoted as with entries

Besides, the shape function is

(6)

It can be proved as Eqs. (4) and (5) that

(7)
(8)
where and are constants independent of h. Equations (7) and (8) indicate that both the determinant and the condition number of are unchanged with respect to h. Consequently, the stabilized IMLS yields more accurate and stable results.

2.3. Error estimate

Along the lines of the error analysis of the IMLS approximation,[36] the following error estimate can be derived for the stabilized IMLS approximation.

3. The IEFG method for linear Robin problem

Since the Robin boundary condition requires fewer restrictions on trial functions, we first consider the following general second-order linear elliptic equation with Robin boundary condition:

(9)
where ς(x) is a symmetric, uniformly elliptic, and bounded matrix, τ(x) ≥ 0 and σ(x) ≥ 0 are given continuous functions, b(x) and are given square-integrable functions, and n(x) is the normal exterior to Γ at the point x. When σ ≡ 0, the Robin problem (9) degenerates to a pure Neumann problem.

The weak solution of the linear Robin problem (9) is to find such that

(10)
where the bilinear form is coercive and continuous on ,
(11)
In order to obtain the meshless numerical solution, let
(12)
where are the stabilized IMLS shape functions given by Eq. (6). Then, the approximation of the variational problem (10) reads to find such that
(13)
which can be discretized into the following linear system
where
(14)
(15)

In Ref. [39], the error estimates of the EFG method are established for the linear Poisson equation. In this paper, the error estimates of the IEFG method will be established for more general linear and nonlinear elliptic equations. As with Theorem 4.1 in Ref. [39], we have the following theorem.

4. The IEFG method for linear Dirichlet problem

Consider the following linear elliptic equation with Dirichlet boundary condition:

(19)
where ς(x) is a symmetric, uniformly elliptic, and bounded matrix, τ(x) ≥ 0 is a given continuous function, and b(x) and are given square-integrable functions.

Due to the use of the weighted least squares technique, the IMLS shape functions lack the delta function property. In this study, the penalty method[1,2] is adopted and then the linear Dirichlet problem (16) is approximated as

which admits the variational problem: Find such that
(20)
where α is a penalty factor, and
(21)

With the stabilized IMLS approximate space given by Eq. (12), the approximation of the variational problem (20) reads to find such that

(22)
which can be discretized into the following linear system
(23)
where u is given in Eq. (14), and
(24)

5. The IEFG method for nonlinear Robin problem

Consider the following nonlinear elliptic equation with Robin boundary condition:

(36)
where is a nonlinear function of the unknown u. Let satisfy the Lipschitz condition with respect to u, i.e.,
(37)
where is a Lipschitz constant.

The weak solution of the nonlinear Robin problem (36) is to find such that

(38)
where the bilinear form is defined in Eq. (11), and

With the stabilized IMLS approximate space given by Eq. (12), the approximation of the variational problem (38) reads to find such that

(39)
which can be discretized into the following nonlinear system
(40)
where u and K are given in Eqs. (14) and (15), respectively, and

Since is a nonlinear term of u, equation (40) is a nonlinear algebraic system and cannot be solved directly. In this research, the system is solved by the subroutine ‘fsolve’ of MATLAB.[43] This algorithm is a variation of Newton’s method, which uses a finite difference approximation to the Jacobian and takes precaution to avoid large step sizes or increasing residuals.[44] In Ref. [44], a boundary type method was developed for nonlinear elliptic problems. The method in Ref. [44] is different from the presented IEFG method. However, the subroutine ‘fsolve’ can be used in the two methods to solve the final nonlinear algebraic systems.

6. The IEFG method for nonlinear Dirichlet problem

Consider the following nonlinear elliptic boundary value problem:

(47)
where is a nonlinear function of the unknown u and satisfies Eq. (37).

By introducing a penalty factor α, we can approximate problem (47) as

which admits the variational problem: Find such that
(48)
where the bilinear form is defined in Eq. (21), and

With the stabilized IMLS approximate space given by Eq. (12), the approximation of the variational problem (48) reads to find such that

(49)
which can be discretized into the following nonlinear system:
where u and are given in Eqs. (14) and (24), respectively, and

7. Numerical examples

In this section, numerical results are presented for several one-dimensional problems. In all examples, the quadratic basis function is employed. Besides, to evaluate the integrations using the standard Gauss quadrature, a uniform cell structure with one node per cell is used, and each cell contains five Gauss points. Moreover, the weight function is chosen as the quartic spline, i.e.,

where and the support size is taken to be 4h.

7.1. Linear Robin problem

Consider the following linear Robin problem:

where and are chosen such that the analytical solution is .

The numerical and analytical values of u and its derivative are presented in Fig. 1 with 33 nodes. We can find that the numerical solutions are in excellent agreement with the analytical solutions.

Fig. 1. (color online) Numerical and analytical solutions of (a) u and (b) for the linear Robin problem.

To study the convergence, figure 2(a) gives the log–log plot of the error with respect to the nodal spacing. Clearly, the error decreases along with the decrease of the nodal spacing, and the IEFG using the stabilized IMLS obtains monotonic convergence. The experimental convergence rate approximates to 2 and agrees well with the theoretical rate. For comparison, the error of the IEFG using the original IMLS is also given in the figure. The comparison shows that the IEFG using the stabilized IMLS has higher computational precision and better numerical stability than the IEFG using the original IMLS. In addition, the curves of the error against the CPU time are shown in Fig. 2(b). It is found that the IEFG using the stabilized IMLS has greater computational efficiency than the IEFG using the original IMLS.

Fig. 2. (color online) Comparison of (a) convergence and (b) computational efficiency of the IEFG using the original IMLS and the stabilized IMLS for the linear Robin problem.
7.2. Linear Dirichlet problem

Consider the following linear Dirichlet problem:

where is chosen such that the analytical solution is .

For investigating the influence of the penalty factor, we display the relationship between the penalty factor α and the error of the IEFG in Fig. 3. In this analysis, four kinds of nodal arrangements of 81, 161, 321, and 641 nodes are used. We can observe that a too small or too big penalty factor increases the numerical error. For a too small penalty factor, Theorem 3 indicates that the error will be dominated by the second right-hand term of Eq. (25) with the theoretical convergence rate . Clearly, the left-hand half of Fig. 3 shows that the experimental convergence rate matches the theoretical one. On the other hand, for a too big penalty factor, Theorem 3 indicates that the error will be dominated by the third right-hand term of Eq. (25) and the error increases as the penalty factor increases. The right-hand half of Fig. 3 also confirms this theoretical result.

Fig. 3. (color online) Convergence of the IEFG with respect to the penalty factor α.

From Fig. 3 we can also find that the optimal value of the penalty factor α is influenced by the nodal spacing h. When the penalty factor is chosen as given in Eq. (35), the log–log plot of the error with respect to h is shown in Fig. 4. To study how the parameter affects the convergence of the IEFG, the error is shown for various . From Fig. 4 we can see that the experimental convergence rate approximates to 2 in all cases and thus the convergence is little influenced by the parameter . These numerical results also match the theoretical results.

Fig. 4. (color online) Convergence of the IEFG with respect to the nodal spacing for variable penalty factors.
7.3. Nonlinear Robin problem

Consider the following nonlinear Robin problem:

where and are chosen such that the analytical solution is .

The numerical results of the IEFG for u and its derivative are shown in Fig. 5 together with the analytical solutions. The results are obtained using 21 nodes on the domain. We can observe from this figure that the numerical solutions are in excellent agreement with the analytical solutions.

Fig. 5. (color online) Numerical and analytical solutions of (a) u and (b) for the nonlinear Robin problem.

The log–log plot of the error with respect to the nodal spacing h is given in Fig. 6. From this figure, we can observe that the experimental convergence rate approximates to 2 and also confirms the theoretical rate.

Fig. 6. (color online) Error and convergence for the nonlinear Robin problem.
7.4. Nonlinear Dirichlet problem

Consider the following nonlinear Dirichlet problem:

where is chosen such that the analytical solution is .

Figure 7 exhibits the analytical and numerical solutions for the potential u and its derivative . In this analysis, the penalty factor is taken as , and 21 nodes are used on the domain. Again, the numerical solutions capture the behavior of the analytical solutions very well.

Fig. 7. (color online) Numerical and analytical solutions of (a) u and (b) for the nonlinear Dirichlet problem.

To demonstrate the convergence of the IEFG for nonlinear Dirichlet problems, we present the log–log plot of errors with respect to the nodal spacing in Fig. 8. It can be seen again that the experimental results match the theoretical bounds.

Fig. 8. (color online) Error and convergence for the nonlinear Dirichlet problem.
8. Conclusion

By discussing the condition number and the determinant of the moment matrix in the IMLS approximation, shifted and scaled orthogonal basis functions are developed to stabilize the IMLS approximation. Compared with the original IMLS approximation, the stabilized IMLS approximation is expected to give more accurate and stable results. Then, the theoretical error of the IEFG method is analyzed for both linear and nonlinear elliptic boundary value problems. For linear and nonlinear Robin problems, the error bound of the numerical solution is proportional to the nodal spacing. For linear and nonlinear Dirichlet problems, since the penalty method is adopted to impose Dirichlet boundary conditions, the error bound is related to the nodal spacing and the penalty factor. Theoretical and numerical analysis indicates that satisfactory results can be achieved by choosing the penalty factor proportional to a negative power of the nodal spacing. Numerical examples confirm the theoretical analysis, and show that the IEFG method using the stabilized IMLS has greater computational efficiency and better numerical stability than the IEFG method using the original IMLS.

Reference
[1] Cheng Y M 2015 Meshless Methods Beijing Science Press in Chinese
[2] Liu G R 2009 Meshfree Methods: Moving Beyond the Finite Element Method Boca Raton CRC Press
[3] Lei J M Peng X Y 2016 Chin. Phys. B 25 020202
[4] Lancaster P Salkauskas K 1981 Math. Comput. 37 141
[5] Wang Y C Li X L 2014 Chin. Phys. B 23 090202
[6] Li X L Li S L 2013 Chin. Phys. B 22 080204
[7] Cheng Y M Li J H 2005 Acta Phys. Sin. 54 4463 in Chinese
[8] Ren H P Cheng Y M Zhang W 2010 Sci. China: Phys. Mech. Astron. 53 758
[9] Cheng Y M Bai F N Liu C Peng M J 2016 Int. J. Comput. Mater. Sci. Eng. 5 1650023
[10] Sun F X Wang J F Cheng Y M Huang A X 2015 Appl. Numer. Math. 98 79
[11] Deng Y J Liu C Peng M J Cheng Y M 2015 Int. J. Appl. Mech. 7 1550017
[12] Cheng Y M Bai F N Peng M J 2014 Appl. Math. Model. 38 5187
[13] Wang J F Sun F X Cheng Y M 2012 Chin. Phys. B 21 090204
[14] Li X L 2015 Appl. Math. Model. 39 3116
[15] Chen S S Wang J Li Q H 2016 Chin. Phys. B 25 040203
[16] Lu Y Y Belytschko T Gu L 1994 Comput. Meth. Appl. Mech. Eng. 113 397
[17] Cheng Y M Peng M J 2005 Sci. China: Phys. Mech. Astron. 48 641
[18] Liew K M Cheng Y M Kitipornchai S 2006 Int. J. Numer. Meth. Eng. 65 1310
[19] Li X L Zhang S G 2016 Appl. Math. Model. 40 2875
[20] Mirzaei D Schaback R Dehghan M 2012 IMA J. Numer. Anal. 32 983
[21] Li X L Li S L 2016 Comput. Math. Appl. 72 1515
[22] Li X L Wang Q Q 2016 Eng. Anal. Bound. Elem. 73 21
[23] Liu W K Jun S Zhang Y F 1995 Int. J. Numer. Meth. Fluids 20 1081
[24] Liu W K Li S Belytschko T 1997 Comput. Meth. Appl. Mech. Eng. 143 113
[25] Belytschko T Krongauz Y Organ D Fleming M Krysl P 1996 Comput. Meth. Appl. Mech. Eng. 139 3
[26] Jin X Z Li G Aluru N R 2001 Comput. Model. Eng. Sci. 2 447
[27] Zhang Z Liew K M Cheng Y M Li Y Y 2008 Eng. Anal. Bound. Elem. 32 241
[28] Zhang Z Liew K M Cheng Y M 2008 Eng. Anal. Bound. Elem. 32 100
[29] Zhang Z Wang J F Cheng Y M Liew K M 2013 Sci. China: Phys. Mech. Astron. 56 1568
[30] Zhang Z Hao S Y Liew K M Cheng Y M 2013 Eng. Anal. Bound. Elem. 37 1576
[31] Zhang L W Deng Y J Liew KM 2014 Eng. Anal. Bound. Elem. 40 181
[32] Peng M J Li R X Cheng Y M 2014 Eng. Anal. Bound. Elem. 40 104
[33] Zhang Z Li D M Cheng Y M Liew K M 2012 Acta Mech. Sin. 28 808
[34] Cheng R J Wei Q 2013 Chin. Phys. B 22 060209
[35] Cheng R J Cheng Y M 2016 Chin. Phys. B 25 020203
[36] Li X L Chen H Wang Y 2015 Appl. Math. Comput. 262 56
[37] Li X L Zhu J L 2009 J. Comput. Appl. Math. 230 314
[38] Li X L 2011 Appl. Numer. Math. 61 1237
[39] Li X L 2016 Appl. Numer. Math. 99 77
[40] Cheng R J Cheng Y M 2008 Acta Phys. Sin. 57 6037 in Chinese
[41] Cheng R J Cheng Y M 2011 Acta Phys. Sin. 60 070206 in Chinese
[42] Li X L Zhang S G Wang Y Chen H 2016 Comput. Math. Appl. 71 1655
[43] Mathews J H Fink K K 2004 Numerical Methods using MATLAB Englewood Cliffs Prentice Hall
[44] Li X L Zhu J L 2009 Eng. Anal. Bound. Elem. 33 322